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A method to calculate NMR J-coupling constants from first principles in extended systems 
is presented. It is based on density functional theory and is formulated within a planewave- 
pseudopotential framework. The all-electron properties are recovered using the projector augmented 
wave approach. The method is validated by comparison with existing quantum chemical calcula- 
tions of solution-state systems and with experimental data. The approach has been applied to verify 
measured J-coupling in a silicophosphate structure, SisO(P04)6. 



I. INTRODUCTION 

Nuclear Magnetic Resonance (NMR) allows informa- 
tion to be relayed through magnetic nuclei in a non- 
destructive and powerful approach to structural elucida- 
tion. It is a fundamental tool in a broad range of scientific 
disciplines and is a cornerstone of modern spectroscopy. 
NMR spectra yield a wealth of information, the most 
commonly reported property being the chemical shift. 
This parameter relates an externally applied magnetic 
field to the resulting change in the local electronic en- 
vironment of the magnetic nuclei, thereby providing key 
insight into the underlying atomic structure. 

NMR J-coupling or indirect nuclear spin-spin coupling 
is an indirect interaction of the nuclear magnetic mo- 
ments mediated by the bonding electrons. It is mani- 
fested as the fine-structure in NMR spectra, providing a 
direct measure of bond-strength and a map of the con- 
nectivities of the system. The J-coupling mechanism is 
an essential component of many NMR experiments.— 

In solution-state, J-coupling measurements can often 
be obtained from one dimensional spectra where the mul- 
tiplet splitting in the peaks is clearly resolved. How- 
ever, in the solid-state this is not the case as these split- 
tings are typically obscured by the broadenings from 
anisotropic interactions. Fortunately this technical chal- 
lenge has not prevented the determination of J-coupling 
in the solid-state, as recent work employing spin-echo 
Magic Angle Spinning (MAS) techniques 2 - has resulted 
in accurate measurements of J-coupling in both inor- 
ganic materials 3 ^ 5 -! 6 - and molecular cry s t als 1 1 1 1 In 
combination with the advances in solid-state experi- 
ments, there has also been an increased interest from the 
biomolecular community as J-coupling has been found 
to be a direct measure of hydrogen bond strength ] 11 ! 12 ! 13 
Both of these factors have provided a strong impetus to 
develop first principles approaches to compute the NMR 
J-coupling constants in order to support experimental 
work, particularly for solid-state systems. 

For finite systems, NMR parameters, including both 
chemical shifts and J-couplings, can be routinely calcu- 



lated using traditional quantum chemistry approaches 
based on localised orbitals. 14 ' 15 Such calculations have 
been widely applied to assign the solution-state NMR 
spectra of molecular systems and establish key confor- 
mational and structural trends^ In particular NMR J- 
couplings have been used to quantify hydrogen bonding^ 
in biological systems. In order to apply these techniques 
to solid-state NMR, it is necessary to devise finite clus- 
ters of atoms which model the local environment around 
a site of interest in the true extended structure. While 
this has led to successful studies of NMR chemical shifts 
in systems such as molecular crystals*** supra-molecular 
assemblies^ and organo-metallic compounds^ it is clear 
that there are advantages in an approach that inherently 
takes account of the long-range electrostatic effects in ex- 
tended systems. 

This observation has led to the recent development 
of the Gauge Including Projector Augmented Wave 
(GIPAW)2i method which enables NMR parameters 
to be calculated at all-electron accuracy within the 
planewave-pseudopotential formalism 2 ^ of density func- 
tional theory (DFT). The technique has been applied, 
in combination with experimental NMR spectroscopy, to 
systems such as minerals ; 23 ' 24 ! 25 glasse s 26 ' 27 and molecu- 
lar crystals^ 2 ^ 

In this paper we introduce a theory to compute 
NMR J-Couplings in extended systems using periodic 
boundary conditions and supercclls with the planewave- 
pseudopotential approach. Like the GIPAW approach to 
calculating NMR chemical shifts, our method is formu- 
lated within the planewave-pseudopotential framework 
using density functional perturbation theory (DFPT). 
We use the projector-augmented-wave 3 ^ (PAW) recon- 
struction technique to calculate J-couplings with all- 
electron accuracy. 

In the following section we discuss the physical mecha- 
nism of the indirect spin-spin interaction, the basis of the 
PAW approach and the supercell technique. In Sections 
IIIII and IIVI we show how the J-coupling tensor maybe be 
calculated using PAW and DFPT. The method has been 
implemented in a parallel plane- wave electronic structure 
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code and we discuss details of the implementation and 
provide validation results in Section |Vl 



II. NMR J-COUPLING 

We consider the interaction of two nuclei, K and L, 
with magnetic moments, /xk and /xl, mediated by the 
electrons. The first complete analysis of this indirect cou- 
pling was provided by Ramse y 32 ^ 33 who decomposed the 
interaction into four distinct physical mechanisms; two 
involving the interaction of the nuclear spins through the 
electron spin and two through the electron charge. In the 
absence of spin-orbit coupling i.e, for relatively light el- 
ements, the charge and spin interactions can be treated 
separately. We can write the magnetic field at atom L 
induced by the magnetic moment of atom K as 
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tl = Rl — r, where Rl is the position of nucleus L, fj,Q 
is the permeability of a vacuum and 5 is the Dirac delta 
function. 

/xk interacts with the electron spin through a mag- 
netic field generated by a Fermi-contact term, which is 
due to the finite probability of the presence of an elec- 
tron at the nucleus, and a spin-dipolar interaction. Both 
of these terms give rise to a first order spin magnetisa- 
tion density, m^^r). This magnetisation density then 
induces a magnetic field at the receiving nucleus by the 
same mechanisms, which in this case are given respec- 
tively by the first and second terms of Eqn. [1] The in- 
teraction between /^k and the electronic charge gives rise 
to an induced current density. To first-order this is given 
by j^( r ) an( i can be divided into a paramagnetic and a 
diamagnetic contribution. 

The J-coupling tensor, J lk> between L and K, can 
be related to the induced field by 



b^RlH 
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where 7l and 7k are the gyromagnetic ratios of nuclei 
L and K. Although the physical interpretation of J- 
coupling is simplified by considering the interaction in 
terms of a responding and a perturbing nucleus, it is 
a symmetric coupling and cither atom L or K can be 
considered as the perturbing site. Experimental interest 
is focused primarily on the isotropic coupling constant, 
Jlkj which is obtained from the trace of J lk and is 
measured in Hz. The superscript, n, denotes the order 
of the coupling in terms of the number of bonds sepa- 
rating the coupled nuclei. In a typical NMR experiment 



J-coupling can be measured across a maximum of three 
bondsi 34 ' 35 In this paper we concentrate solely on obtain- 
ing the isotropic or scalar value. 

To calculate J kl we obtain m' 1 ) and jW within den- 
sity functional perturbation theory using a planewave ex- 
pansion for the wavefunctions with periodic boundary 
conditions and pseudopotentials to represent the ionic 
cores. The use of pseudopotentials generates a compli- 
cation as the J-coupling tensor depends critically on the 
wavefunction in the regions close to the perturbing and 
receiving nuclei, precisely the regions where the pscudo- 
wavefunctions have a non-physical form. To compensate 
for this we perform an all-electron reconstruction of the 
valence wavefunctions in the core region using BlochPs 
projector augmented wave (PAW) schemed 

Within this scheme, the expectation value of an oper- 
ator O, applied to the all-electron wavefunctions, is 
expressed in terms of the pseudised wavefunctions, 
as: (ip\0\ip) = (il>\0\ip). Here, for an all-electron lo- 
cal or semi-local operator, O, the corresponding pseudo- 
operator, O, is given by 



(1) 
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R labels the atomic site, or augmentation region, and 
n and m are composite indexes which account for the 
angular momentum channels and the number of projec- 
tors. |</>R,ri) are the all-electron partial waves obtained as 
eigenstates of an atomic calculation within r c , the pseu- 
dopotcntial core radius and \4>B.,n) are the correspond- 
ing pseudo partial waves. |pR, n ) are the localised pro- 
jectors which weight the superposition of partial waves 
where (pk, re |0R/, m ) = S RR 'S nm . The PAW method has 
been used to calculate several all-electron properties from 
pscudopotcntial calculations including: EPR hypcrfinc 
parameters^ electric field gradient tensors 3 -! and Elec- 
tron Energy Loss Spectroscopy^ 

To calculate J-couplings in the solid-state using peri- 
odic boundary conditions, the perturbing nucleus can be 
viewed similar to a defect in a defect calculation. This 
allows us to use the standard technique of constructing 
supercells from the unit cell which are large enough to in- 
hibit the interaction between the periodic defects or per- 
turbations. This corresponds to extending the system- 
size to facilitate the decay of the induced magnetisation 
and current densities within the simulation cell. Figure. 
[T]is a schematic of a 2 x 2 x 2 superccll constructed from 
eight unit cells. The perturbing atom now lies at the cor- 
ner of a much larger cell which decreases the interaction 
between the perturbation and its periodic image. This 
approach works very well for localised properties such as 
J-coupling. To calculate the J-coupling for molecules, we 
use a vacuum superccll technique. In both cases, the J- 
couplings must be converged with respect to the cell-size. 
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FIG. 1: Schematic of the supercell technique. The unit cell is 
on the left, a 2x2x2 supercell of the unit cell is on the right. 
This supercell doubles the distance between the perturbing 
atom (black) and its periodic image in the next cell. 

III. THE SPIN MAGNETISATION DENSITY 

We now obtain the contribution to the J-coupling ten- 
sor which arises from the interaction of the nuclear spins 
mediated by the electron spin. We first obtain an ex- 
pression for the pseudo-Hamiltonian in the presence of a 
perturbing nuclear spin and show how it can be used to 
obtain the induced magnetisation density. We then use 
the magnetisation density to calculate the magnetic field 
induced at the receiving nucleus. 

A. Pseudo-Hamiltonian 

The all-electron Hamiltonian for a system containing 
N magnetic moments which interact through the elec- 
tron spin, S, is expanded to first order in the magnetic 
moment of the perturbing site, /^k, to give 

H = ip 2 + V (0) (r) + V (1) (r) + H SD + H FC . (4) 

where 

H SD =.g/JS-B| D , (5) 

and 

H FC = g(3S ■ Bl c . (6) 

B^ D is the magnetic field generated by a spin-dipole in- 
teraction, 

tdSD _ Mo 3r K (r K • Mk) - r^ K l 

and B^ c is the Fermi contact interaction, 

B£ C = y^K)/^. (8) 

We have defined rx = Rk — r, where Rk is the position 
of nucleus K and X is the identity matrix. Here V^ ^ (r) is 
the ground-state all-electron local potential and V*- 1 -* (r) is 
the corresponding first order variation. The latter term 



is due to the change in magnetisation density induced 
by [j>k- The perturbation does not give rise to a first 
order change in the charge density and so there is no 
change in the Hartree potential for linear order. This 
can be understood by considering the effect of time re- 
versal; the charge density is even under time inversion, 
while the spin-magnetisation and magnetic field change 
sign. As a result the perturbation does not induce a first 
order change in the charge density, however there is a first 
order induced spin-magnetisation density. V^^r) there- 
fore accounts for the first-order variation of the exchange- 
correlation term which we label as K^J . 

We now use the PAW transformation (Eqn. [3]) to ob- 
tain the pseudo-Hamiltonian. As Eqn.[3]does not contain 
any field dependence it is sufficient to apply it to each 
term of Eqn. [4] individually. The pseudo-Hamiltonian at 
zcroth order in the perturbation is 

H(°) = ip 2 +V loc (r) + ^VS, (9) 

R 

where Vi oc (r) is the local part of the pseudopotentials 
which includes the self-consistent part of the Hamilto- 
nian and is the non-local part which is given by 
V S = E„, m \PR,n)a* m (pR, n \- a* m are the strengths 
of the nonlocal potential in each channel at each ionic 
site. Collecting terms to linear order in the perturbation 
gives, 

H« =H« +H SD +H FC . (10) 

HgD describes the spin-dipolar interaction induced by /xk 
and is given by 

H SD = g(1S ■ B| D + g/3S ■ ABjp. (11) 

The first term on the right-hand side is the all-electron 
operator and the second term is the augmentation to this, 

AB| D = IPR."> I (^«|B| D |^R, m ) (12) 

-(0 R ,„|B| D |0 Rim )](p Rim |, 

with R = Ra'- In Eqn. [TT] we have only included the 
augmentation of the spin-dipolar operator at the site of 
the perturbing atom. This on-site approximation is fully 
justified given the localised nature of this operator. 

Hpc is the Fermi-contact operator and can be con- 
structed in a similar manner to the spin-dipole operator 
giving an all-electron and an augmentation contribution. 
However, as the Fermi-contact operator contains a Dirac 
delta-function and is therefore localised within the aug- 
mentation region, Hpc can be simplified considerably. 
The pseudo-partial waves and projectors, ^ n (j> n \<p n ), 
form a complete set which enables us to rewrite the all- 
electron contribution in terms of the pseudo-partial waves 
within the augmentation region and so we can equiva- 
lently express the operator as 

H FC = g/3S ■ |PR,n)(<W|B£ c |0R, ra >(pR, m |, (13) 

n, m 
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where R = Rk- This form is more suitable for a practical 
calculation as it avoids an explicit representation of the 
delta-function. 



B. Magnetisation Density 

To construct the magnetisation density, we define 
m^(r) to be the linear response to the magnetic field, 
Bi induced along the direction by the spin-dipolar 
and Fermi-contact interactions. The total magnetisation 

density is obtained as m' 1 ' = Ei=a; y z m i 1 ^( r ): the sum 
over the cartesian directions. By choosing as the spin 
quantisation axis, is diagonal in the spin-up and 

spin-down basis. The eigenstates of H~(°) + H^ are also 
eigenstates of u, • S and so the magnetisation density is 
parallel to giving; 



ra.f'(r)]ui 



(r)uj 



(14) 



Here 



mV (r) = 5 /3 [n$ (r) - (r)] - 2g(3n^ (r), (15) 

where g and (5 were defined previously and n^l is the 
induced density for either spin-up or spin-down. The 
simplification of the magnetisation density in this way 
is a consequence of time reversal symmetry, namely the 
absence of a first-order charge density. This means that 
the spin-up and spin-down ground-state wavefunctions 



are equivalent, so that 



l« } >- 



Also, the linear 



variation of the wavefunctions induced by the spin mag- 
netisation are related through \tpig) 
Within PAW, Eqn. [T5l becomes 



4 3 /3Re]>>J 1 ) |r><r|< ) } 
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Re signifies taking the real component, ) are the 

eigenstates of the unperturbed Hamiltonian, H' ', {ip^J) 
are the perturbed pseudowavefunctions and o indexes the 
occupied bands. The first term on the right hand side of 
Eqn. [in] is the pseudo-magnetisation density fn- 1 , and 
the second term is the corresponding augmentation. For 
simplicity, we drop the spin indexing on the ground-state 
wavefunctions from now on as the spin-dependence enters 
only through the perturbation. 

To calculate \tp^g) we employ a Green's function 
method where 



I ; (0) \ / / (0) I 

WJ) = Q(e)\^) = V K ^ W irt*)- (17) 



is the first order Hamiltonian given by Eqn. [TOl with 
the spin quantised along the direction. G(e) is the 
Green's function, e G and e e are the eigenvalues of the oc- 
cupied and empty bands. Rather than explicitly sum over 
the empty states, we project onto the occupied bands by 
multiplying Eqn. [T7] through by (e — H^ ) ) . We define 

7> = E e l$ 0) ><$ 0) | = 1 - E l^ 0) )(^ 0) | and rewrite 
Eqn. [T7] as 



(e -H(°))|^)=PH«|^°)). 



(18) 



This is then solved using a conjugate gradient minimisa- 
tion scheme, with an additional self-consistency condition 
to account for the dependence of Hi* on the spin-density. 
For a more detailed account of this type of approach, see 
Ref. M. 



C. Induced Magnetic Field 

The induced magnetic field at atom L, and subse- 
quently the J-coupling between L and K, due to the spin 
magnetisation is obtained by combining Eqns. [Tl and [TBI 
to give 



B% (R L ) = B^(R L ) + AB^(R L ) + AB^(R L ), 

(19) 

where we have taken advantage of the linearity of Eqn. [JJ 
to yield three separate terms. The first term is the pseudo 
spin-dipolar contribution and the second term is the spin- 
dipole augmentation. The notation used here implicitly 
assumes that a rotation over the spin-axis has been per- 
formed. 

In practise, the pseudo-spin dipole term can be con- 
structed from the Fourier transform of Eqn. [JJ 
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(20) 

where m^^G) is the Fourier transform of the pseudo- 
magnetisation density and G is a reciprocal space lattice 
vector. The G = term is neglected as the cell size 
is large compared with the strength of the perturbation 
which is small. 

The induced magnetic field at atom L is then recov- 
ered by performing a slow inverse Fourier transform at 
the position of each responding nucleus. The spin dipole 
augmentation term is obtained by using an on-site ap- 
proximation and evaluating terms of the form, 

AB«(R L ) = .9/?£RcE (V^IABE !^). (21) 



ABjp is defined in Eqn. [JJ] but now the subscript indi- 
cates the responding rather than the perturbing nucleus. 
To evaluate this term and Eqn.ll2| we note that \<p n ) can 
be decomposed into the product of a radial (|R n ;)) and an 
angular (|Yz m )) term. BfP can also be rewritten as the 
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product of a radial and angular component such that the 
computation of the augmentation term involves the on- 



site calculation of (Rni\-^t\Rn'i')(Yr, 



T\Y m ). The 



latter quantity reduces to an integral over spherical har- 
monics given by the Gaunt coefficients. 

The Fermi-contact contribution, ABp^ , is obtained by 
following the same argument used in constructing Hpc 
(Eqn. I13[) and is given by 



AB«(R L ) 



(22) 



which is the all-electron reconstruction of the induced 
magnetisation density at the responding nucleus. 

IV. CURRENT DENSITY 

We now obtain the contribution to the J-coupling ten- 
sor arising from the interaction of the nuclear spins me- 
diated by the electron charge current. The derivation of 
the current density is similar to that of the magnetisation 
density and much of the notation is conserved through- 
out. We first obtain an expression for the pseudo- 
Hamiltonian in the presence of a perturbing nuclear spin 
and show how it can be used to obtain the induced cur- 
rent density. We then use this current density to calculate 
the magnetic field induced at the receiving nucleus. 

A. Pseudo-Hamiltonian 

The all-electron Hamiltonian for a system of N mag- 
netic nuclei which interact through the nuclear vector 
potential, A(r) = £ N A N (r) = MY^/™, C£m be 
expanded to first order in /xk to give, 



Hi. 



f 



p 2 +V ( °)(r)+H, 



where 



H AK = £p.A K (r), 



(23) 



(24) 



and V^ ) (r) is the ground-state local potential. The per- 
turbation does not induce a first order change in ci- 
ther the charge or magnetisation densities and so un- 
like Eqn. |H there is no linear variation to the self- 
consistent potential in Eqn. [23l We have used the sym- 
metric gauge for the vector potential and taken the nat- 
ural choice of gauge-origin; namely that for the Nth nu- 
clear spin the gauge origin is the Nth atomic site giv- 
ing A N (r) = l/2B(r) x r N . This gauge-choice preserves 
the translational invariancc of the system and is much 
simpler than in the otherwise analogous case of NMR 



chemical shielding. In the latter situation, due to the use 
of finite basis sets a rigid translation of the system in 
the uniform external magnetic field introduces an addi- 
tional phase factor. For the planewave-pseudopotential 
approach the problem was solved by Pickard and Mauri 
with the development of the Gauge-Including Projec- 
tor Augmented- Wave (GIPAW) approach, an extension 
which is unnecessary here. 

To obtain the pseudo-Hamiltonian we now apply the 
PAW transformation of Eqn. [3] to Eqn. [23J The zeroth 
order term is again given by Eqn. [S] and the first order 
term by 



Mo r K x p 
4tt MK ' |r K l 3 



Mo 

4-7T 



n.rn 



fR,n I 



FX 



;I<?W 



(25) 



The first term on the right-hand side is the all-electron 
component and the second term is the augmentation 
which is constructed at the site of the perturbing nu- 
cleus. Lk = i"k x p is the angular momentum operator 
centred on the perturbing atomic site. 



B. Current Density 

The current density operator, J(r), is given by the sum 
of a paramagnetic and a diamagnctic term, 

J(r) = J p (r)+J d (r), (26) 

where the paramagnetic term is given by 

JP(r) = -[p|r)(r| + |r)(r|p]/2, (27) 

and the diamagnetic term is 

J d (r) = -A(r)|r)(r|. (28) 

If we consider the current due only to the perturbing 
nucleus, and with our atomic choice of gauge origin, the 
diamagnetic term can be written as 



J|(r) 



Mo Mk x r K 

47r rt 



r)(r| 



(29) 



By applying Eqn.[3jto both Eqns. [27l and [29l we obtain 
the pseudo-current density operator within PAW 



J(r) = JP(r) + 4(r) + [AJP(r) + AJ d (r)] , 

where the paramagnetic augmentation operator is 

AJP(r) = \PR, n )[(<fa.,n\J p \<fci, m ) 



(30) 



R,,n,m 

-(^R,n|J P |0R,m) ] (pk,i 



(31) 
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and the corresponding diamagnctic operator is 
AJ^(r) = J2 [(MJi\<hi,m) 

R, , n , in 

-(^R,n|JKl?R,m) ] (PR,m|- (32) 

Arranging terms in J(r) to zcroth and linear order in /xk 
gives 



and 



j(°)(r) = JP(r) + AJ p (r) 



jW(r)=j|(r)+A4(r) 



(33) 



(34) 



Using Eqns. [33] and [34] we are now able to obtain the 
first-order induced current density, which within density 
functional perturbation theory is given by 

j«(r) = 2^2Re(^ )|j(o)|^(i) ) + (^0)^(1)^(0)^ 

O 

= jW(r)+jW(r). (35) 

Here is the unperturbed wavefunction, {ipo 1 ^} is 

the perturbed wavefunction and o indexes the occupied 
bands. The first term on the right hand side is the para- 
magnetic contribution to the induced current and the sec- 
ond term is the diamagnetic contribution. The first order 
wavefunction is again obtained using Eqn. 1181 where IT 1 ' 
is now given by Eqn. 1251 



C. Induced Magnetic Field 

We can now combine Eqns. Q] and [35] to calculate the 
J-coupling between nuclei L and K arising from the mag- 
netic field induced by the orbital current. The magnetic 
field at nucleus L due to the orbital current can be ex- 
pressed as the sum of 4 terms, 

B« (R L ) = B« (R L ) + > (R L ) + AB« (R L ) 

+AB^(R L ), (36) 

the pseudised contributions from the paramagnetic and 
diamagnetic currents and their respective augmentation 
terms. 

To calculate the pseudised contributions to the current 
density we obtain the Fourier transform of b£, (Rl) and 
B^Rl) giving 



B (1) (G) 



iGxj« d (G) 



G 2 



(37) 



where jpV d (G) is the Fourier transform of either the para- 
magnetic or diamagnetic current. To obtain the induced 



field at the atom site L we perform a slow Fourier trans- 
form of Eqn. [37] We again note that the G=0 contri- 
bution toB( 1 )(G) is neglected as the contribution is ex- 
pected to be small. The augmentation to the paramag- 
netic current is calculated using an on-site approximation 
(R = R L ) with 



47T 



bn,n\—!r\4>R,r 



f>R,n\^-\^R,m) 



(PR, m |^ 1} ), (38) 



where L is the angular momentum operator evaluated 
with respect to the augmentation regions. The augmen- 
tation to the diamagnetic current is given by 



AB 
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This is much more difficult to evaluate than any other 
term as an on-site approximation cannot be used due to 
the presence of the rx position vector within the augmen- 
tation summation. However, previous quantum chemical 
studies have shown that the overall contribution to the 
J-coupling from the diamagnetic term is very small com- 
pared with any of the other three contributions. In light 
of this, we have neglected this term in our current imple- 
mentation. 



V. RESULTS 

We have implemented our theory into a parallelised 
plane-wave electronic structure code4£ The ground- 
state wavefunctions and Hamiltonian arc obtained self- 
consistently after which the isotropic J-coupling constant 
is calculated using the outlined approach. In our im- 
plementation we use norm-conserving Troullier-Martins 
pscudopotcntials4i For the all-electron reconstruction we 
used two projectors per angular momentum channel. 

In the following sections we compare our approach 
to existing quantum chemistry approaches which use lo- 
calised basis sets and with experiment. 



A. Molecules 

To validate our method we have calculated isotropic 
coupling constants for a range of small molecules^ and 
compared them with experiment. There are several 
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FIG. 2: J-couplings calculated for a set of molecules. Both 
the MCSCF and BLYP results were taken from Ref.[H All of 
the experimental values were also taken from this paper. The 
PW(PBE) results are from the present work. The lines are 
obtained from a linear regression of the calculated values with 
experiment. All values are quoted in 10 19 T 2 J _1 , the unit of 
the reduced coupling constant. 



studies of calculated J-couplings for small molecules re- 
ported in the literature, using a variety of theoretical ap- 
proaches, see Ref. [HI and references therein. We compare 
our results to calculations presented by Lantto et ai^ 3 . 
which were obtained within DFT using the BLYP func- 
tional and with the Multi-configurational Self-Consistent 
Field (MCSCF) approach. For consistency we use the 
molecular geometries reported in their work. 

To obtain the isotropic J-coupling we use a supercell 
of size 1728 A 3 for each molecule with the exception of 
benzene which required a larger cell-size of 3375 A 3 . The 
exchange-correlation was approximated by GGA-PBE 44 
and an energy cut-off of 80 Ry was imposed on the 
planewave expansion. All calculations sample the Bril- 
louin zone at the gamma-point and used norm-conserving 
Trouillier-Martins pseudopotentials.— The calculated J- 
coupling against experiment for the molecules are shown 
in Fig. [5] alongside the results of a linear regression for 
each set of data. These results are presented as reduced 
spin-coupling constants, which are given by Klk = 
2 5 J LK . and so are independent of nuclear species. The 

R7L7K ' 1 1 

graph indicates an excellent overall agreement with both 
experiment and the other approaches. The accuracy of 
the planewave PBE calculations is comparable with the 
all-electron BLYP results, with regression coefficients of 
0.92 and 1.08 respectively. The correlation coefficients 
for the BLYP, MCSCF and PW data are 0.97, 0.97 and 
0.99 respectively, suggesting a smaller random error in 
the planewave approach. Unsurprisingly, the most ac- 
curate couplings are given by MCSCF, which provides 
a more comprehensive description of electron correlation 
but is computationally more demanding than DFT. 

In Table U we present the J-coupling values calcu- 
lated for benzene. The results compare favourably with 



TABLE I: J-coupling[Hz] in benzene. PW(PBE) labels the 
current planewave approach. The BLYP, MCSCF and ex- 
perimental values were taken from Ref. |43j. D, P, FC and 
SD label the diamagnetic, paramagnetic, Fermi-contact and 
spin-dipolar contributions respectively. All values are in Hz. 



Jkt, Method 



D 



J cc 



•I, 1 



Jibe- PW(PBE) 
BLYP 
MCSCF 
J cc PW(PBE) 
BLYP 
MCSCF 
PW(PBE) 
BLYP 
MCSCF 
PW(PBE) 
BLYP 
MCSCF 
J CH PW(PBE) 
BLYP 
MCSCF 
J CH PW(PBE) 
BLYP 
MCSCF 
J CH PW(PBE) 
BLYP 
MCSCF 



FC SD Total Expt[Hz] 



0.2 
0.2 
0.4 
0.0 
0.0 
-0.2 
0.0 
0.0 
0.1 
0.5 



-0.1 



-6.9 
-6.8 
-6.6 
0.0 
0.1 
0.0 
0.6 
0.5 
0.4 
0.9 



58.2 
63.8 
75.1 
-1.2 
0.0 
-3.7 
7.3 
8.4 
16.8 
132.5 



0.2 5.1 



-0.2 -1.1 6.0 



-0.1 



0.3 -0.4 



1.9 
1.0 
1.5 
-0.3 
-0.5 
-1.1 
1.1 
1.5 
1.8 
-0.2 



0.1 



0.0 



0.0 



53.4 
58.2 
70.9 
-1.5 
-0.4 
-5.0 
9.0 
10.4 
19.1 
133.7 
155.2 
185.1 
5.3 
1.1 
-9.8 
4.7 
7.4 
12.9 
-0.2 
-1.3 
-6.1 



55.8 



-2.5 



10.1 



158.3 



1.0 



7.6 



-1.2 



both the existing approaches and with experiment. The 
MCSCF approach systematically overestimates the J- 
coupling for both Jch and Jcc compared with experi- 
ment. This is due to the use of a restricted basis set 
which was necessary given the size of the system, for fur- 
ther details see Ref. |45[ 

The decomposition of the J-coupling into the four com- 
ponents serves as an illustration of the relative strengths 
of each contribution and the trends over several bonds. 
Lantto et al have only presented this separation for Jcc • 
It is clear that the Fermi-contact is the dominant mech- 
anism in the coupling and that the diamagnetic compo- 
nent is consistently the smallest and is often negligible. 



B. Crystals 

Due to the difficulties encountered measuring J- 
coupling in solid-state systems there are very few val- 
ues to be found in the literature that are suitable for 
validation of our approach. Recently Coelho et alA pro- 
vided an estimate for the two bond coupling between 
29 Si and 31 P pairs in the silicophosphatc SisO^O^e- 
Subsequently this was followed by a more accurate 
determination^ which identified the four Si-O-P cou- 
plings. We have calculated NMR chemical shifts and 
J|_ _ Si for SisO(P04)6 to validate our approach. 

The structure of Si 5 0(P0 4 ) 6 is trigonal (a=7.869A, 
c=24.138A, 36 atoms per primitive cell) and contains 
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TABLE II: Calculated J-coupling for silicophosphate TABLE III: Calculated NMR chemical shifts^ and J-coupling 
Si 5 0(P0 4 ) 6 using the unit cell and two supercells constructed for silicophosphate Si 5 0(P0 4 )6- The experimental values are 
with 2x1x1 and 2x2x1 unit cells. in brackets and were taken from Ref. [j. 



Coupling 


lxlxl 


2x1x1 


2x2x1 


Coupling 


31 P [ppm] 29 Si [ppm] 


Calc. [Hz] 


^Jp-Oa-Sii 


-17.37 


-17.07 


-17.12 


Jp-Oa-Si! 


-47.4 (-43.8) -214.8 (-213.3) 


-17.12 (15±2) 


Jp-o 2 -si 2 


-16.16 


-16.18 


-16.26 


T 2 

Jp_0 2 -Si 2 


-218.7 (-217.0) 


-16.26 (14±2) 


2 T 

JP-0 6 -Si 2 


-1.30 


-1.20 


-1.17 


T 2 

Jp_0 5 -Si 2 


-218.7 (-217.0) 


-1.17 (4±2) 


2 T 

JP-0 4 -Si 3 


-13.83 


-14.18 


-14.13 


Jp_0 4 -Si3 


-128.6 (-119.1) 


-14.13 (12±2) 



one unique P site and three incquivalent Si sites. Two of 
these Si sites are 6-fold coordinated, Sii and Si2, and the 
third site, S13, is four- fold coordinated. Sii is bonded to 
six equivalent oxygen atoms, Si2 is bonded to six oxygen 
atoms which are comprised of two distinct sites. Si3 is 
bonded to three equivalent oxygen atoms and one oxy- 
gen from an S13-O tetrahedron. Thus there is one 31 P 
chemical shift, three 29 Si chemical shifts and four unique 
2 J P _o-Si couplings. 

We obtained the structure from the Chemical Database 
Service at Daresbury.— Prior to calculating the NMR 
parameters, we performed a full geometry optimisa- 
tion on the structure, using a planewave cut-off of 70 
Ryd and norm-conserving pseudopotentials. The GGA- 
PBE^i exchange-correlation functional was used and a 
Monkhorst-Pack k-point grid with a maximum of 0.1 
A" 1 between sampling points. We calculated the NMR 
chemical shifts using the GIPAW— approach with the 
same parameters used for the geometry optimisation. 
The J-coupling between P and Si was obtained using 
the approach outlined above. A slightly higher maxi- 
mum planewave energy (80Ryd) was required to give J- 
couplings converged to within 0.1Hz. 

We tested the convergence of the induced magnetiza- 
tion density and current density with respect to super- 
cell size. The results of these calculations for three cell 
sizes are presented in table [TTJ From this we can see 
that both the induced magnetization and current densi- 
ties have decayed substantially within the single unit cell. 
The largest of these calculations (144 atoms) was paral- 
lelised over 16 dual-core AMD processors and took 45 
hours to run. The groundstate calculation took approxi- 
mately 14 hours and the J-coupling terms; Fermi-contact, 
spin-dipolar and orbital, required 3.5 hours, 15.5 hours 
and 11.4 hours respectively. 

The results for the 2x2x1 cells are presented in com- 
parison with experiment in Table IIIIl From Table IIIII 
it is clear that the calculated J-couplings are in excel- 
lent agreement with experiment and fully reproduce the 
surprisingly large spread in the J-coupling values. Our 
calculations verify the novel experimental work and also 
identify the sign of the couplings which are not deter- 
mined by the experimental spin-echo based approaches. 

The NMR chemical shifts are also in good agreement 
with experiment, particularly for 29 Si. For both 29 Si 
and 31 P the difference between the calculated and ex- 
perimental values is a very small fraction of the total 



TABLE IV: Decomposition for the J-coupling in Si50(P0 4 )e. 
D is the diamagnetic term, P is the paramagnetic term, SD 
is the spin-dipolar and FC is the Fermi-contact. 



Coupling 


D 


P 


SD 


FC 


Total 


2j P _o 3 -Sii 


-0.05 


-0.27 


-0.03 


-16.77 


-17.12 


2 T 

J P _0 2 -Si 2 


-0.02 


-0.50 


-0.23 


-15.51 


-16.26 


2 T 

Jp-o 5 -si 2 


-0.10 


-0.07 


0.18 


-1.18 


-1.17 


2 T 

JP-0 4 -Si 3 


-0.09 


-0.49 


0.23 


-13.79 


-14.13 



shift range. We note that our assignment of the three Si 
sites in Si50(P04)6 agrees with the assignment based on 
experimental intensities as discussed by Coelho et a&. In 
Table ITVl we present the decomposition of the silicophos- 
phate J-coupling into their constituent terms. As with 
benzene, the Fermi-contact is found to be consistently 
the largest component while the diamagnetic and spin- 
dipolar contributions are very small. 

VI. CONCLUSIONS 

We have developed an all-electron approach for calcu- 
lating NMR J-coupling constants using plancwavcs and 
pseudopotentials within DFT. Our method is applicable 
to both solution and solid state systems using supercell 
techniques. We have validated our theory against ex- 
isting quantum chemical approaches and experiment for 
molecules. We have calculated the J-coupling between Si 
and P in a silicophosphate polymorph, for which we have 
determined the sign of the coupling. 

Given the recent experimental interest in J-coupling, 
we expect that our approach will prove useful in deter- 
mining both the range and strength of coupling in sys- 
tems not yet investigated and whether or not such cou- 
plings can feasibly be determined by experiment. By 
combining J-coupling calculations with computations of 
other NMR parameters, there now exists a comprehen- 
sive set of computational tools to complement experimen- 
tal understanding and design. 
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